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We perform axisymmetric simulations of tlie magnetorotational collapse of very massive stars in 
full general relativity. Our simulations are applicable to the collapse of supermassive stars with 
masses M > 10^ M© and to very massive Population III stars. We model our initial configurations 
by n = 3 polytropes, uniformly rotating near the mass-shedding limit and at the onset of radial 
instability to collapse. The ratio of magnetic to rotational kinetic energy in these configurations 
is chosen to be small (1% and 10%). We find that such magnetic fields do not affect the initial 
collapse significantly. The core collapses to a black hole, after which black hole excision is employed 
to continue the evolution long enough for the hole to reach a quasi-stationary state. We find that 
the black hole mass is Mh = 0.95M and its spin parameter is Jh/M^ — 0.7, with the remaining 
matter forming a torus around the black hole. The subsequent evolution of the torus depends on 
the strength of the magnetic field. We freeze the spacetime metric ("Cowling approximation") and 
continue to follow the evolution of the torus after the black hole has relaxed to quasi-stationary 
equilibrium. In the absence of magnetic fields, the torus settles down following ejection of a small 
amount of matter due to shock heating. When magnetic fields are present, the field lines gradually 
coUimate along the hole's rotation axis. MHD shocks and the magnetorotational instability (MRI) 
generate MHD turbulence in the torus and stochastic accretion onto the central black hole. When the 
magnetic field is strong, a wind is generated in the torus, and the torus undergoes radial oscillations 
that drive episodic accretion onto the hole. These oscillations produce long-wavelength gravitational 
waves potentially detectable by the Laser Interferometer Space Antenna (LISA). The final state of 
the magnetorotational collapse always consists of a central black hole surrounded by a coUimated 
magnetic field and a hot, thick accretion torus. This system is a viable candidate for the central 
engine of a long-soft gamma-ray burst. 

PACS numbers: 04.25.Dm,97.20.Wt,97.60.-s 



I. INTRODUCTION 

Population III stars born with zero metallicity com- 
prise the first generation of stars. It is believed that 
their formation causes the reionization of the universe 
and terminates the "dark ages" (see, e.g., Ill] and ref- 
erences therein). The disruption of Pop III stars fol- 
lowing nuclear burning may be responsible for the small 
metal abundance observed in later generations of stars 
(e.g. Pop II stars). Simulations of the collapse of pri- 
mordial molecular clouds suggest that Pop III stars tend 
to be massive. Masses in the range between lOOM© and 
IOOOMq are not uncommon Some of these calcu- 
lations suggest that the initial mass function for Pop III 
stars has a bimodal distribution, with peaks at ~ lOOM© 
and 1-2Mq |3i]- Stars with masses between 14OM0 and 
260Mq encounter a pair instability and are likely to be 
completely disrupted by nuclear-powered explosions 4]. 
The recent observation of the peculiar Type Iln super- 
nova SN2006gy in NGC1260 points to the possibility that 
such a disruption can occur in a massive star (> I2OM0) 
even during the current epoch Q . For stars with masses 
above 260M©, the explosive nuclear burning is unable to 
reverse the implosion and the stars are likely to collapse 
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directly to black holes f3| . 

Growing evidence indicates that supermassive black 
holes (SMBHs) with masses in the range 10*^ - IO^^Mq 
exist and are the engines that power active galactic nuclei 
(AGNs) and quasars [1, There is also ample evidence 
that SMBHs reside at the centers of many, and perhaps 
most, galaxies ^ , including the Milky Way @ . The high- 
est redshift of a quasar discovered to date is zqso = 6.43, 
corresponding to QSO SDSS 1148-1-5251 Accord- 
ingly, if they are the energy sources in quasars (QSOs), 
the first SMBHs must have formed prior to zqso — 6.43, 
or within t = 0.87 Gyr after the Big Bang in the concor- 
dance ACDM cosmological model. This requirement sets 
a significant constraint on black hole seed formation and 
growth mechanisms in the early universe. Once formed, 
black holes grow by a combination of mergers and gas 
accretion. 

The more massive the initial seed, the less time is re- 
quired for it to grow to SMBH scale and the easier it 
is to have a SMBH in place by z > 6.43. One possible 
progenitor that readily produces a SMBH is a supermas- 
sive star (SMS) with M > lO^Mg [1, SMSs can 
form when gaseous structures build up sufficient radia- 
tion pressure to inhibit fragmentation and prevent nor- 
mal star formation; plausible cosmological scenarios have 
been proposed that can lead to this situation [l^. Al- 
ternatively, the seed black holes that later grow to be- 
come SMBHs ma y o riginate from the collapse of Pop III 
stars < lO^Af0 [1^. To achieve the required growth 
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to ~ lO^M© by ZQso ^ 6.43, it may be necessary for 
gas accretion, if restricted by the Eddington limiting lu- 
minosity, to occur at low efficiency of rest-mass to ra- 
diation conversion (< 0.2). Recent relativistic simula- 
tions [ii, M, M 

show that accretion onto a rotating 
black hole that has reached spin equilibrium does occur at 
low efficiency in a magnetized disk with turbulence driven 
by the magnetorotational instability (MRI) 17, 18, f9]. 
Such accretion may enable a Pop 111 seed to achieve the 
necessary growth by z = 6.43 [20] ■ But it may be more 
difhcult to use the Pop 111 seeds to explain the origin of 
the first generation of the SMBHs should quasars be de- 
tected at redshifts significantly higher than zqso = 6.43. 

Recent simulations of binary black hole mergers sug- 
gest that gravitational radiation reaction can induce a 
large kick velocity (> f 000 km/s) in the remnants follow- 
ing mergers [1^. These large kick velocities may pose a 
great hazard for the growth of black hole seeds to SMBHs 
by z ^ 6 [iSl , but such large kicks are possible only if the 
spins of the black hole binary companions are appreciable 
and their masses are comparable. Determining the spins 
of the seed black holes formed from collapse and track- 
ing their subsequent evolution via accretion and minor 
mergers [l^ [25!| is therefore important for estimating the 
kick velocities following major mergers. 

It is very likely that massive Pop III and SMSs are 
rotating and have magnetic fields. A SMS does not 
reach sufficiently high temperature for nuclear burning 
to become important before the onset of the general 
relativistic radial instability (26j . Quasistatic contrac- 
tion driven by radiative cooling will spin up the star to 
the mass-shedding limit, provided that viscosity and/or 
magnetic fields are sufficient to maintain uniform ro- 
tation [2^. The star will then evolve secularly along 
the mass-shedding limit, simultaneously emitting electro- 
magnetic radiation, matter, and angular momentum (see, 
e-g- [2l|,|2^[23|). After reaching the onset of radial insta- 
bility, the star collapses on a dynamical timescale. Dur- 
ing the collapse, the rotation becomes differential and the 
rotational and magnetic energies are both amplified. The 
black hole that forms will be rotating and surrounded by 
a magnetized accretion disk. A qualitatively similar final 
fate should characterize a massive Pop III star > 260Mq. 

Shibata and Shapiro performed the first full general 
relativistic (GR) simulation of the collapse of a very mas- 
sive, rotating star [28|. They modeled the massive star 
as a uniformly rotating n = 3 polytrope spinning at the 
mass-shedding limit at the onset of radial instability. A 
massive star is supported largely by thermal radiation 
pressure and is adequately modeled by an n = 3 poly- 
trope. They terminated their simulation soon after the 
black hole formed because of numerical inaccuracies as- 
sociated with the spacetime singularity that inevitably 
forms inside the black hole. They estimated the final 
state of the system using semi-analytic methods (see 
also d^llOl). They concluded that, independent of the 
initial mass M of the progenitor star, the mass of the 
black hole that forms is Mh ~ 0.9M and the hole spin 



parameter is Jh/M'^ ~ 0.75. The remaining gas forms a 
rotating torus around the nascent black hole. 

In this paper, we first repeat the full GR (axisymmet- 
ric) simulation performed by Shibata and Shapiro of mas- 
sive star collapse to the appearance of a black hole. We 
then employ the technique of black hole excision [si], [s^l 
to continue the evolution. We are able to follow the 
spacetime evolution for another 200M by this means. By 
this time, the central black hole and the spacetime metric 
have both settled down to a quasi-stationary state. We 
find that the mass and spin parameter of the final black 
hole are Ah w 0.95M and Jh/M^^ « 0.7. These results 
are close to the semi- analytic estimates in [2^, [2^, [30| . 
The torus surrounding the black hole continues to evolve 
long after the black hole has settled down. In order to 
study the subsequent evolution of the torus, we adopt the 
"Cowling approximation" whereby we freeze the metric 
at i ~ 150Af after the excision and continue to evolve the 
system for another 2000 Af. We find that a small amount 
of material I0~"^Af ) is ejected from the system due to 
shock heating, and the torus relaxes to a dynamical equi- 
librium state ~ lOOOA/ after the formation of the central 
black hole. 

Next, to study the important role of magnetic fields, 
we add a small, seed poloidal magnetic field to the initial 
rotating star and follow the collapse once again. We con- 
sider two different strengths of the seed magnetic fields 
(models SI and S2). The initial magnetic energy Ai is 
1% of the initial rotational kinetic energy T for model 
SI, and 10% of T for model S2. Since T/\W\ = 0.009, we 
have A^/|W^| ^ 1 in both models, where W is the gravi- 
tational potential energy. Hence in both cases the mag- 
netic fields represent small perturbations to the dynam- 
ics of the initial star. During the collapse, the frozcn-in 
poloidal field is amplified as a result of compression. The 
development of differential rotation generates a toroidal 
field due to magnetic winding. However, we find that 
magnetic fields do not affect the collapse significantly be- 
fore the formation of the central black hole. The final 
mass and spin parameter of the black hole are about the 
same as in the unmagnetized case. But magnetic fields 
do affect the evolution of the torus significantly. Mag- 
netic fields intensify the outflow of the ejected material. 
The outflow also lasts longer than in the unmagnetized 
case. As the torus evolves, magnetic fields are collimated 
along the black hole's rotation axis. For model SI, MHD 
shocks and the MRI in the torus create turbulence, which 
leads to stochastic accretion of material from the torus 
to the central black hole. For model S2, a strong wind 
is generated (possibly by the magneto-centrifugal mech- 
anism (Hi) during the period ~ 900Af-1200Af follow- 
ing central black hole formation. This wind induces a 
radial oscillation of the torus, which leads to episodic 
accretion of material to the central black hole, and long- 
wavelength gravitational radiation potentially detectable 
by the Laser Interferometer Space Antenna (LISA). 

The final state of the magnetorotational collapse con- 
sists of a central black hole surrounded by a collimated 
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magnetic field and a massive, hot, accretion torus. These 
features provide the essential ingredients for generating 
ultrarelativistic jets at large distance. Our simple equa- 
tion of state (EOS) is a reasonable approximation for the 
collapse of SMSs, but our omission of neutrino emission 
and other microphysics is certainly not adequate to cap- 
ture all of the physical processes occurring during the 
collapse of massive Pop III and Pop I/II stars. Neverthe- 
less, we expect that the black hole-torus remnant that we 
find will be qualitatively similar to the remnants formed 
from these progenitors if they are rotating rapidly at the 
onset of collapse. The reason is that these stars may 
also be crudely modeled by n « 3 polytropes initially 
and their EOSs may also be represented by an adiabatic 
F « 4/3 law during collapse (see [30l|). 

Our simulations may also help explain the formation of 
the central engine in the coUapsar model of long-soft 
gamma-ray bursts (GRBs). In addition, some GRBs ob- 
served at very high redshift might be related to the gravi- 
tational collapse of very massive Pop III stars [1^ . Hence 
our simulations may also provide insights into the forma- 
tion of GRB central engines arising from these stars. 

The reminder of this paper is organized as follows. In 
Secini we briefly describe the mathematical formulation 
of the Einstein-Maxwell-MHD coupled equations and nu- 
merical techniques used to solve them. We then describe 
our initial data and computational setup in Sec. IIIII We 
present our numerical results in Sec. IIVI and provide a 
summary of our simulations in Sec. |Vl Throughout this 
paper, we adopt geometrical units in which G = 1 = c, 
where G and c denote the gravitational constant and 
speed of light, respectively. Cartesian coordinates are 
denoted by x'^ — {x,y,z). The coordinates are ori- 
ented so that the rotation axis is along the z-direction. 
We define the coordinate radius r = + + z'^, 

cylindrical radius uj — \/ x"^ +2/^, and azimuthal angle 
ip — ta.n~^{y/x). Coordinate time is denoted by t. Greek 
indices fi,!/, - ■ ■ denote spacetime components (t, x, y, z), 
small Latin indices denote spatial components 

(x, y, and z). 



II. FORMULATION 

A. Basic equations and numerical metliods 

The formulation and numerical scheme for our 
GRMHD simulations are the same as those reported 
in 36|, to which the reader may refer for details. Here 
we briefly summarize the method and introduce our no- 
tation. 

We use the 3-1-1 formulation of general relativity and 
decompose the metric into the following form: 

ds^ = -a^dt^ + -f,j{dx' + P'dt){dx' + f3^dt) . (1) 

The fundamental variables for the metric evolution 
are the spatial three-metric 7ij and extrinsic curva- 



ture Kij. We adopt the Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) formalism [33] to evolve jij and Kij . 
In this formalism, the evolution variables are the confor- 
mal exponent = ln7/12, the conformal 3-metric 7ij = 
e~^'^7ij, three auxiliary functions P' = — 7*-' j , the trace 
of the extrinsic curvature K, and the tracefree part of the 
conformal extrinsic curvature Aij = e~'^'^{Kij — ^ijK /Z). 
Here, 7 = det(7ij). The full spacetime metric g^i, is re- 
lated to the three- metric 7^1^ by 7^1/ = Qfiiy +Tii_in^, where 
the future-directed, timelike unit vector n'^ normal to the 
time slice can be written in terms of the lapse a and shift 
/3* as = a-^{l,-f3^). 

The Einstein equations are solved in Cartesian coor- 
dinates. In this paper, we assume both equatorial and 
axisymmetry so we only evolve the region with a; > 
and z > 0. We adopt the Cartoon method [ll] to impose 
axisymmetry in the metric evolution, and use a cylindri- 
cal grid to evolve the MHD and Maxwell equations. As 
for the gauge conditions, we adopt the hyperbolic driver 
conditions as in [s^] to evolve the lapse a and shift /3'. 

The fundamental variables in ideal MHD are the rest- 
mass density po, specific internal energy e, pressure P, 
four- velocity u^, and magnetic field B^^ measured by a 
normal observer moving with a 4-velocity (note that 
B'^Ufj^ = 0). The ideal MHD condition is written as 
U/^iF'^" = 0, where -F'"^ is the electromagnetic tensor. 
The tensor F'"' and its dual in the ideal MHD approxi- 
mation are given by 
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(2) 
(3) 



where e^^a/^ is the Levi-Civita tensor. Here we have in- 
troduced an auxiliary mag netic 4- vector &^ = B'^^-^/V^, 
where is the magnetic field measured by an observer 
comoving with the fluid and is related to i?^ by 



(ti) 



{d''^+u''u^)B'' 



(4) 



The energy-momentum tensor is written as 



where TjJ^"^ and denote the fluid and electromag- 
netic pieces of the stress-energy tensor. They are given 
by 



(6) 



(7) 



where h = 1 + e -I- P/po is the specific enthalpy, and 
&^ = 6^6^,. Hence, the total stress-energy tensor becomes 



Tf_,y = [poh + b )Uf^Uy + ( f + y ) - bf_,by . (8) 
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In our numerical implementation of the GRMHD and 
magnetic induction equations, we evolve the densitized 
density p*, densitized momentum density Si, densitized 
energy density f, and densitized magnetic field B^. They 
are defined as 



(9) 
(10) 

(11) 
(12) 



During the evolution, wc also need the three- velocity = 

The GRMHD and induction equations are written in 
conservative form for variables p*. Si, f, and and 
evolved using a high-resolution shock-capturing (HRSC) 
scheme. Specifically, we use the monotonized central 
(MC) scheme [s^ for data reconstruction and the HLL 
(Harten, Lax and van-Leer) scheme [40l | to compute 
the flux. The magnetic field has to satisfy the no 
monopole constraint diB"^ = 0. We adopt the flux- 
interpolated constrained transport (flux-CT) scheme [i^ 
to impose this constraint. This scheme guarantees that 
no magnetic monopoles will be created in the computa- 
tional grid during numerical evolution. At each timestep, 
the primitive variables (po? v^) must be computed from 
the evolution variables (p*, f , St). This is done by numer- 
ically solving the algebraic equations (|9|)- (fTT|l together 
with an EOS P = P{po,e). 

As in many hydrodynamic simulations in astrophysics, 
we add a tenuous "atmosphere" that covers the compu- 
tational grid outside the star. The atmospheric rest-mass 
density is set to « 10~^°pc(0) before the black hole forms, 
where Pc(0) is the initial rest-mass central density of the 
star. In the excision evolution where the system consists 
of a central black hole and a surrounding torus, the max- 
imum density in the torus is ~ lOOpc(O), and we set the 
atmosphere density to 10~^pc(0). 

The codes used here have been tested in multiple rel- 
ativistic MHD simulations, including MHD shocks, non- 
linear MHD wave propagation, magnetized Bondi accre- 
tion, and MHD waves induced by linear gravitational 
waves [36i. We have also compared this code with the 
GRMHD code developed independently by Shibata and 
Sekiguchi by performing simulations of the evolution 
of magnetized, differentially rotating, relativistic, hyper- 
massive neutron stars lisl. l46l , and of magnetorotational 
collapse of stellar cores (47f7^We obtain good agreement 
between these two independent codes. 



B. Equation of state 



where pressure is dominated by thermal radiation. For a 
Pop I/II star, which has smaller mass, the pressure of the 
pre-collapse core is dominated by the relativistic degen- 
erate electron pressure, which is also well-approximated 
by a r = 4/3 EOS. During the collapse, the EOS stiff- 
ens when the density exceeds nuclear density pnuc ~ 
2 X 10^^ g cm~^. However, if the mass of the collaps- 
ing core exceeds a critical value Merit, the black hole 
horizon appears before the star reaches the nuclear den- 
sity. In this case, the stiffening of the EOS has no effect 
on the collapse. To estimate Merit, consider the collapse 
of a uniform density dust sphere (Oppenheimer-Snyder 
collapse). A horizon appears when the areal radius of 
the sphere reaches R = 2M. At this time, the density 
is po = 3Af/(47ri?3) 1.7 x 1016(Mq/M)2 g cm'^. Set- 
ting po = pnuc, gives Merit ~ lOM©. For a SMS, the 
mass is much larger than Merit. For a Pop III star of 
mass M = 300Mq, the mass of the collapsing core is 
180Mq 3, which is still much larger than Merit. Hence 
the r = 4/3 EOS is also a good approximation during the 
entire collapse phase for very massive Pop III stars. For 
Pop I/II stars, on the other hand, the core mass is less 
than 2Mq and a more realistic EOS is required in the late 
stages. In addition, neutrino emission and transport are 
also important to the dynamics of the collapse for these 
stars. Neutrino generation and transport also play a role 
in the collapse of Pop HI stars ^, i5Q] , but are probably 
not dynamically important for the most massive progen- 
itors or for SMSs because of their low temperature and 
density. 



C. Diagnostics 

During the evolution, we monitor the L2 norm of the 
Hamiltonian and momentum constraints as in [4^ . We 
find the violation of the constraints is at most a few per- 
cent before excision. After the excision, the constraints 
can rise to w 10%. We terminate the excision evolution 
before the constraints reach ^20%. 

We also compute the rest mass A/q, ADM mass M 
and angular momentum J during the evolution. They 
are computed by the following volume integrals: 



Mo = 
M = 

J = 



167r ^ 247r ' 



1 piffcp 
1671 



1 



167r 



(13) 

(14) 
(15) 



In this paper, we adopt the simple n = 3 (F = 4/3) 
polytropic EOS to construct the initial model and P — 
(F - l)poe (F-law EOS) during the evolution. This EOS 
is a good approximation for the pre-collapse core of a 
massive Pop HI star [41,] or the bulk of a SMS [HI, [Hi, 



where F^fe is the Christoffel symbol and R is the Riemann 
scalar associated with ^ij. Note that the above formula 
for J is only valid in an axisymmetric spacetime jsij . The 
rest mass Mg is conserved as a result of the baryon num- 
ber conservation. Angular momentum J is conserved in 
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axisymmetry, as gravitational radiation carries no angu- 
lar momentum. However, M is not conserved since grav- 
itational radiation carries energy and propagates off the 
computational grid. We find that M remains constant to 
within 2%. Our finite difference scheme guarantees that 
Mo and J computed from the above volume integrals are 
conserved to machine precision provided that no fluid 
leaves the computational grid. However, we perform sev- 
eral regriddings during the calculation (see Section III) 
and these leave behind a few percent of Mq and J in the 
outermost layers. 

During the excision evolution, we compute the rest 
mass Mdisk and angular momentum Jdisk of the disk out- 
side the black hole by computing integrals (fT5|) and (fT5|) 
over the volume outside the apparent horizon. The irre- 
d ucible m ass Mi„ of the black hole is given by Mi„ — 

A/l&TT, where A is the surface area of the apparent 
horizon. Since J is conserved, we can compute the black 
hole's angular momentum Jh by 



above expressions reduce to 



Jh ~ J — Jlo 



Jd 



isk 



(16) 



where Jioss is the loss of angular momentum as a result 
of regriddings and matter leaving the grid. The black 
hole's mass is then computed from the formula 



(17) 



which is exact for a Kerr spacetime, and is in accord 
with the formula derived using the isolated and dynami- 
cal horizon formalism [5^. 

At At K, 150A/ after the excision evolution, wc find 
that the spacetime becomes nearly stationary. In this 
case, the energy E is approximately conserved thereafter, 
where 



E 



j "V7 



(18) 



We can then define the fluxes of rest mass, energy, and 
angular momentum across any closed two-dimensional 
surface S" in a time slice: 



(19) 
(20) 
(21) 



FE{r) = - j aT\(fi:, , 
Fj{r) ^ j^aT^Sll, , 



J r— const 

FE{r) = -(f dAa^T\ 

J r— const 



F.Ar) 



dAa^T^^, 



(23) 
(24) 



(25) 



where dA = sm9d9d(t). The total energy flux Fe is 
very close to the rest-mass flux Fm since Fe is primar- 
ily composed of the rest-mass energy flow. Thus, we 
define another energy flux by subtracting the rest-mass 
flow: Fe ~ Fe — Fm- We note that Fe contains kinetic, 
thermal, electromagnetic, and gravitational potential en- 
ergy fluxes. If i^e > at sufficiently large radius, an 
unbound outflow (overcoming gravitational binding en- 
ergy) is present. 

Another method to determine whether a fluid particle 
is unbound is to compute ut- In a stationary spacetime, 
the value of Ut of a particle moving on a geodesic is con- 
served. If the particle is unbound, the radial velocity 
> and —ut — l/\/l — v'^ > 1 at infinity. Hence 

and Ut are useful diagnostics to determine if the fluid 
element is unbound, provided that the fluid motion is 
predominantly ballistic and pressure and electromagnetic 
forces can be neglected. This is usually the case in the 
low-density region. 

During the excision evolution, the MRI may develop in 
the torus surrounding the central black hole. The growth 
time (e-folding time) and wavelength of the fastest- 
growing MRI mode can be roughly estimated by the fol- 
lowing formulae derived in linear perturbation theory in 
Newtonian gravitation [l^ \^ : 



^MRI 



2|9r2/91nro 
2ttv 



1- 



- 

\2n 



-1/2 



(26) 
(27) 



where is angular velocity, = /y^AirpQ is the z- 
component of the Alfven speed and 



1 diw'^n^'^ 



dz 



1/2 



(28) 



is the epicyclic frequency. 



m. INITIAL DATA AND GRID SETUP 
A. Initial Data 



where 



d Yii — -iijkdx^ A dx 



(22) 



and Cijk = n^e^ijk is the Levi-Civita tensor associated 
with the three-metric 7^ . If S" is a sphere of radius r, the 



We model the pre-collapse star as a uniformly rotating 
star satisfying the n = 3 polytropic EOS P = Kp^^ . We 
set K = 1 in our code. As explained in [ill, our result 
can be scaled to arbitrary values of K or, equivalently, the 
ADM mass M; only nondimensional ratios are invariant. 
For example, M cx AT^/^, B cx AT"^/^, po oc K~^,... etc. 
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We use the same initial model as in [28*1, whereby 
the star is rotating near the mass-shedding limit with 
r/|VF| = 0.009. The equatorial radius of the star is 
= 640Af 10'^(M/10^Mo)km. This is the con- 
figuration where the polytrope is on the verge of ra- 
dial instability against gravitational collapse due to gen- 
eral relativity [20|. The central density of the star is 
Pc = 1O3(M/1O4M0)-2 g cm-3. For a SMS with mass 
M > IO^Mq, this general relativistic instability triggers 
the collapse, as opposed to microphysical processes such 
as pair instability. 

We add a small amount of poloidal seed magnetic field 
to this equilibrium star, employing a magnetic vector po- 
tential of the form 



I 1/6 

max(pQ 



(29) 



where pcut — 10"'^ /Oc, and Ab is a constant that deter- 
mines the strength of the initial magnetic field. The mag- 
netic field is computed by the formula = t^^^djA^. A 
similar form of the initial magnetic field has been used in 
the study of magnetized accretion disks around a station- 
ary black hole H . IZsj , the collapse of hypermassive neu- 
tron stars [islkGji and the collapse of a magnetized, stel- 
lar core of a high mass star f47j . We choose two nonzero 
values of the constants A^ so that the values of initial 
magnetic energy 



(30) 



are 1% and 10% of the initial rotational kinetic energy 
(corresponding to A^/|W^| = 9 x 10^^ and 9 x lO^'^). 
Hence adding this seed magnetic field causes only a slight 
perturbation to the star. We label these two models as 
SI and S2, respectively. We also study the unmagne- 
tized case {A4 = 0, model SO) to compare with the pre- 
vious result reported in [2^ . We can also characterize the 
strength of the magnetic field by the volume-averaged ra- 
tio of gas pressure to the magnetic pressure. Specifically, 
we define /3 = {P)/{Pmag), where the magnetic pressure 
is Pmag — and 



id) 



JqdV 



v.. 



Here dV — ^/^d^x is the proper volume element and 
Vs — /p>o volume of the star. This definition 

of [3 is used in in the study of magnetized accre- 
tion disks around central black holes. The value of (3 for 
models SI and S2 are 3700 and 370, respectively. We 
also define the averaged strength of magnetic field B by 
B = ^SttAI /Vs- In cgs units, we find 



B = 3 X 10 



8 ( 



M 



\ lO^Mp 



G 



for model SI and 
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FIG. 1: Initial density contour curves (thick, black) and mag- 
netic field lines (thin, green). The density contour curves 
are drawn for po = 10~'~'''^pc with j — 0, 1, • • • , 10, and 
the poloidal magnetic field lines (for models SI and S2 only), 
which coincide with contours of in axisymmetry, are for 



Ain Air 



c(ji/20) with j — 1^2, - ■ ■ ,19 where pc and A^^ 



denote the central density and maximum value of Aip, respec- 
tively. Note that although the magnitudes of the magnetic 
fields are different for models SI and S2, the field lines have 
the same profile when normalized as described. 



for model S2. 

The strength of the magnetic field inside a Pop III star 
is unknown and is currently not addressed by theoreti- 
cal models dealing with their cosmological formation [i^ . 
Our goal is to determine what effects, if any, magnetic 
fields may have on the eventual collapse of the stars to 
black holes. After all, they are possible progenitors of 
GRBs, and many GRB models require a magnetized disk 
around a black hole. Here, we choose the strengths of the 
seed magnetic field to be sufficiently large for us to per- 
form reliable simulations with limited computational re- 
sources and still be able to resolve the wavelength of the 
fastest growing MRI mode. We note that these magnetic 
field strengths are still quite small dynamically (small 
and large /?), and so the magnetic field does not 
affect the dynamics of the collapse (see Sec. IIVI) . How- 
ever, the post-collapse evolution does depend on the cho- 
sen strengths. In Sec IIVI we will discuss how our results 
may change for even smaller field strengths. 

Following [ill, we induce collapse by depleting 1% of 
the pressure (i.e., P 0.99P) everywhere inside the star. 
The parameters of our models are summarized in TablelH 
The density and magnetic field profiles of our pre-coUapse 
model are shown in Figure [T] 
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TABLE I: Model parameters 



Model 


MjT 


/3 


B X (M/lO*Af0) 


SO 











SI 


0.01 


3700 


3 X 10*G 


S2 


0.10 


370 


lO^G 



B. Grid Setup 

We perform simulations using a cell-centered uniform 
grid with size iV x 3 x iV in x-y-z, covering a computa- 
tional domain A/2 < a; < L- A/2, A/2 <z< L- A/2, 
and —A < y < A. Here, N and L are constants and 
A = L/N. The variables in the y = ±A planes are com- 
puted from the quantities in the y — plane by imposing 
axisymmetry. Since the characteristic radius of the star 
decreases by a factor of ^ 1000 during the collapse (from 
~ 600M to ~ M), using a fixed uniform grid with suffi- 
cient resolution for the entire collapse phase is computa- 
tionally prohibitive. In order to save computational re- 
sources and at the same time ensure adequate resolution 
throughout the simulation, we adopt a regridding tech- 
nique similar to the algorithm described in 2^i|. When 
gravity is weak (in the Newtonian regime), the charac- 
teristic radius of the star is proportional to 1/(1 — ac), 
where ac is the central lapse. We thus use a regridding 
algorithm based on the values of ac- During the early 
stages, the collapse proceeds in a homologous manner. 
We set N = 400 and L = 929M when ac > 0.984. Keep- 
ing N fixed, we decreases L as the collapse proceeds: 
L = 656M when 0.976 < < 0.984, L = 459M when 
0.905 < ac < 0.976. After this stage, the collapse in 
the core is faster than in the outer layers. We increase 
the grid number N and decrease A as follows: TV = 900 
and L = 158M when 0.7 < ac < 0.905, N = 1400 and 
L = 135M when 0.3 < ac < 0.7. In the last stage, the 
star collapses to a black hole. In order to allocate our grid 
more effectively in this last stage, we interpolated the 
data onto a multiple-transition fisheye coordinates (ssj 
when ac < 0.3. 

The multiple transition fisheye coordinates are re- 
lated to the original coordinates through the following 
transformation: 




10 15 20 
t/(l03M) 



30 



FIG. 2: Evolution of central rest-mass density (upper panel) 
and lapse (lower panel) for models SO (black solid lines), SI 
(red dotted lines) and S2 (blue dashed lines). The central 
density is normalized by its initial value Pc(0). Note that 
the results for SO and SI are very close and their lines almost 
overlap. The plots terminate soon after the apparent horizons 
appear. 



118M. When transformed back to the original coordinate 
system, the outer boundary is « 60Af and the resolutions 
are 



A 



0.025M r<4Af 

0.05M 4M < r < llAf 

O.IM llAf <r< 22M 

0.2Af r > 22M 



(34) 



We find that the total rest mass and angular momen- 
tum that are discarded as a result of the regriddings are 
about 1% and 5-8% of their initial values for the models 
considered. 



IV. RESULTS 



x^ = —r{r), 
r 

n 

r{f) = a„f -I- ^ Ki In 

i=l ' 

_ (oj-i - ai)si 
2tanh(foi/si) ' 



cosh[(r + roi)/s,] 
cosh[(f- f(ii)/si] 



(31) 
(32) 

(33) 



where r — ■\/arM-yM-?, f ~ \/ xP' -\- y'^ + ^ n, ai, 
foi and Si are constant parameters. In the last stage of 
collapse (ac < 0.3), we use a cell-centered uniform grid 
with N — 600 in fisheye coordinates with parameters n = 
3, (oo,ai,a2,a3) = (0.125,0.25,0.5,1), (roi, ro2, rosl = 
{31.5M, 59.5M, 81AM), si = S2 = S3 = 5.69Af, and L = 



Figure [5] shows the evolution of central density and 
lapse for the three models. Figures [3]-[5] show the density 
contours and velocity vectors during pre-excision evo- 
lution for models SO, SI and S2 respectively. Poloidal 
magnetic field lines are also shown for models SI and S2 
in Figs. [4] and [5) We see that magnetic fields slightly 
slow down the collapse. As mentioned in Sec. IIIIl the 
collapse proceeds in a homologous manner at the begin- 
ning. When the central lapse decreases to ac < 0.9, 
the central region collapses faster than the outer layers. 
The apparent horizon appears at t — 28280Af for model 
SO, t = 28360Af for SI and t = 29149Af for S2. With- 
out excision, the code becomes inaccurate soon after the 
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FIG. 3: Density contour curves and velocity vector fields in the meridional plane for model SO (pre-excision) . The density 
levels are drawn for po = /5scailO"°-^^ (j = 0-12), where pscai = llpc(O) at t = 25260M, pscai = 340pc(0) at t = 27930, and 
Pscai = lOOOpc(O) at t = 28284M. The thick (red) line near the lower left corner in the far right graph denotes the apparent 
horizon. Note that the scale is different for each time slice to show the central region in detail. 



t = 25420M 



= 28000M 



t = 28364M 




30 40 50 
X/M 



FIG. 4: Density contour curves and velocity vector fields (upper graphs), and magnetic field lines (lower graphs) in the 
meridional plane for model SI. The thick (red) lines near the lower left corner in the far right graphs denote the apparent 
horizon. The density levels are drawn for po = Pscai 

IQ-o-Sj = 0-12), where p^cai = llpc(O) at t = 25420M, p.cai = 340pc(0) 
ai t = 28000, and pscai = lOOOpc(O) at t = 28364M. The poloidal magnetic field lines are drawn as contours of A^, with levels 
given by A^p — (j/20)A^3,max with j — 1-19, where ^,p,max is the maximum value of A^ at the given time. Note that the scale 
is different for each time slice to show the central region in detail. 
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FIG. 5: Same as Fig. H but for model S2. 



formation of the apparent horizon because of the grid 
stretching. 

To continue the evolution, we excise a spherical region 
inside the apparent horizon. We start the excision evo- 
lution a few Ai M after the apparent horizon forms. 
We are able to follow the evolution reliably for another 
~ 200M, after which the Hamiltonian and momentum 
constraints increase substantially. This eventual break- 
down is probably because the metric inside the horizon, 
which is not computed accurately, slowly leaks out to the 
region outside due to superluminal gauge modes. We are 
currently investigating other gauge conditions, as well as 
other techniques to overcome the numerical difRculty. As 
the collapse proceeds, the mass and angular momentum 
of the central black hole increase before settling down 
to quasi-stationary values. Figures [SHI] show the post- 
excision evolution of the black hole's irreducible mass 
Miir, mass Alh, spin parameter Jh/M^, and the rest mass 
of the material outside the apparent horizon Mdisk, for 
the three models. We see that after At ^ 150M, the 
black hole settles down to a quasi-stationary state, with 
Mh « 0.95M and Jh/M^ « 0.7 for aU the three mod- 
els. This result agrees roughly with our earlier simula- 
tions and analytic estimates for unmagnetized collapse 
{Mh « 0.9M and Jh/M^ « 0.75) in [28, 29, 30]. The re- 
maining material, having too much angular momentum, 
forms a torus surrounding the black hole (see Figs. [S]- 
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FIG. 6: Post-excision evolution of the mass Mh, spin parame- 
ter Jh /Mh , and the irreducible mass Min of the central black 
hole, and the rest mass of the disk Afdisk outside the apparent 
horizon for model SO. Time is measured from the beginning 
of excision [U^ = 28284M). 
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FIG. 7: Same as Fig. [6] but for model SI. Excision starts at 
tcx = 28364M. 
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FIG. 8: Same as Fig. [6] but for model S2. Excision starts at 
tex = 29150M. 



[TT|) . Even though the central black hole has settled down 
after ~ 150M, the torus continues to evolve as material 
from the outer layers gradually reaches the central region. 
The dynamical timescale at radius r is tdyn ~ l-K^r^ jM. 
Hence idyn ~ lOOOM at r = 30M, and tdyn ~ 2000M at 
r = 50M. Since the torus extends beyond 50M, we need 
to follow the evolution for at least 2000M. To study the 
subsequent evolution, we adopt the Cowling approxima- 
tion by freezing the metric at t — tex ~ 150M, where iox 
is the time when excision starts. This is a fairly good ap- 



proximation since the material outside the horizon con- 
tributes only ~ 5% of the total mass and so the metric 
is dominated by the central black hole, which has settled 
down. We have compared the results of our Cowling (sta- 
tionary metric) and non-Cowling (dynamic metric) runs 
during the transition interval 150M ^ t ~ tcx 200M 
and find good agreement. 

Figure [5] shows snapshots of density contours and ve- 
locity fields for model SO in the post-excision evolution. 
We terminate the simulation at i — tex = 2200M, where 
most of the dynamical processes have ended. We find 
that an outflow develops at i — tox ^ 120M near the 
horizon and becomes prominent at t — t^x ~ 170M. 
The outflow is due to the fact that material from the 
outer layers arrives into the inner region with a sub- 
stantial amount of angular momentum. When it reaches 
the inner region, the centrifugal barrier prevents it from 
falling into the black hole. The fluid particles move with 
"zoom- whirl" -like trajectories [541] and accumulate near 
the black hole. As more fluid particles arrive and smash 
into interior layers, the fluid heats up and forms a shock, 
which propagates outward and creates an outflow along 
the surface of the torus. While this outflow is still ex- 
panding, we flnd that a secondary, weaker outflow forms 
at t — tox ~ 160Af, which can be seen in the second and 
third plots in Fig. [9l A few more episodes of smaller 
outflow develop as more material from the outer layers 
arrive. However, the process damps hy t > 250M by 
which time most of the material has reached the cen- 
tral region and the residual infalling fluid does not have 
enough momentum to push on the torus and generate 
further outflow. To determine if the outflowing material 
is unbound, we calculate the quantity —Ut- As discussed 
in Sec. [Tll any unbound fluid particle moving in a low 
density region (in which pressure and electromagnetic 
forces are negligible) has —Ut > 1. We flnd that the 
outflow material is indeed unbound, but the total rest 
mass of the unbound fluid is only 10~'^M. The outflow 
reaches the outer boundary of our grid (r « 60M) after 
t — tex ^ 500M. Most of the unbound material leaves the 
grid after t — tex ^ 700M. During this same time, the 
infalling material in the outer region of the torus close 
to the equatorial plane also rebounds outward because of 
the centrifugal barrier. We have checked that this out- 
ward moving fluid remains bound (— Ut < 1), but about 
0.02M of rest mass leaves the grid by the time we termi- 
nate our simulation at t — tex = 2200M. This outward 
moving fluid has too much angular momentum to be able 
to remain in the inner region. The torus in the inner re- 
gion with r < 30M settles down to quasi-equilibrium by 
t > 500M. 

Figures [TUl and [Tll show snapshots of density contours, 
velocity fields and poloidal magnetic field lines for mod- 
els SI and S2 in the post-excision evolution. We find 
similar outflow as in the case SO, but the outflow in S2 
develops at time t — t^x ~ 75M, much earlier than that 
in SO. Unlike SO, the outflow in SI and S2 is generated 
continuously rather than intermittently. The outflow is 
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FIG. 9; Snapshots of density contour curves and velocity vectors in the post-excision evolution for model SO. The contours are 
drawn for po ~ 100/ac(0)10~'''^-' (j = 0-10). The red line denotes the apparent horizon. 



stronger than SO and about 4 x 10~'^M of the rest mass 
becomes unbound for SI and 8.7 x Vd~^M for S2, much 
larger than the case of SO. In the presence of magnetic 
fields, the outflow carries the frozen magnetic field and 
travels outward along the torus's surface. This causes the 
field lines near the boundary of the outflow and torus to 
bend (see Figs. [TUl and [TT|) . This bending amplifies the 
magnetic field in that region and hence the outflow is 
intensified by the extra magnetic pressure. A magnetic 
shock is also generated in the region, which leads to tur- 
bulence in the torus. The bending is more significant in 
SI than in S2. This is because in S2, the magnetic field 
is strong enough to quickly counteract the bending and 
drives more fluid outward. Figure [T^ shows the rest-mass 
flux Fm, energy flux and angular momentum flux Fj 
through a spherical surface of radius 50Af for the three 
models. We see that the outflow is significantly stronger 
in the presence of magnetic fields. Figure [T^ also indi- 
cates that a sustained fiux is present in the time period 
900M-1200M for model S2. We find that this flux is 
not due to outflow generated near the black hole, but 
due to a wind that arises in the middle of the torus. We 
find that in the wind the fluid moves along the mag- 
netic field lines. The inclination angle between the field 
lines and z-axis is between 20° and 40°. This suggests 
that the wind is driven by the magneto-centrifugal mech- 



anism [33|. The outflows in models SI and S2 cause the 
fleld lines to coUimatc along the rotation axis of the black 
hole. For model S2, the outflow and the subsequent wind 
carry away a substantial amount of magnetic energy from 
the torus. At t — ^ox > 1500M, the wind subsides and 
the interior of the remaining torus has a weak magnetic 
fleld. The outflow and wind in model S2 are so strong 
that they perturb the equilibrium of the inner torus and 
causes it to oscillate radially. As in the case of SO, there 
is bound fluid moving out of the grid as a result of cen- 
trifugal bounce in models SI and S2. By the end of the 
simulation if, - t^^ = 2000M), only 0.04M of the rest 
mass remains in the inner torus in model SI and 0.02M 
remains in model S2. 

Figure fT3l shows the rest- mass flux through the appar- 
ent horizon for the three models. For model SO, the in- 
ward flux decreases with time as the torus settles down to 
dynamical equilibrium. Without magnetic fields or vis- 
cosity, there is no dissipation to drive further accretion. 
For model SI, we see that at late time {t — tex ?J 800M) 
material from the torus accretes into the central black 
hole in a stochastic manner. Stochastic accretion is often 
seen in simulations of magnetized accretion disks around 
stationary black holes (see e.g. [H, HI]). This suggests 
that the accretion is due to magnetic-induced turbulence 
in the torus. The turbulence is generated initially by the 
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FIG. 10: Snapshots of density contour curves and velocity vectors (first and third rows), and poloidal magnetic field lines 
(second and fourth rows) in the post-excision evolution of model SI. The contours are drawn for po = 100pc(0)10-° =^^ (j = 0- 
10). The thick (red) line near the lower left corner denotes the apparent horizon. The poloidal magnetic field lines are drawn 
for = (j/20)A^,max with j = 1-19. 




FIG. 11: Same as Fig.[TO]but for model S2. 
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FIG. 12: Rest- mass flux Fm, energy flux and angular mo- 
mentum flux Fj through a spherical surface of radius 50Af 
for models SO (black solid lines), SI (red dotted lines) and S2 
(blue dashed lines). 
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FIG. 13: Rest-mass flux Fm through the apparent horizon 
for models SO (black solid lines), SI (red dotted lines) and S2 
(blue dashed lines). 



magnetic shock as a result of the outflow, and is then 
sustained by the MRI. To verify that we are able to re- 
solve the MRI, we compute the wavelength of the fastest 
growing MRI mode Amri using Eq. (|27|) . We find that 
Amri/A >15 in some region near the equatorial plane, 
where A is our grid spacing. This suggests that the MRI 
can be resolved in our simulation. For model S2, the 



radial oscillation of the inner torus causes episodic accre- 
tion into the central black hole. When the torus swings 
away from the black hole, no accretion occurs. Accretion 
resumes when the torus swings towards the black hole. 
This explains the episodic mass accretion pattern seen in 
Fig. [T31 The small accretion rate in the figure is due to 
accretion from the atmosphere. 

When the magnetic field strength is much smaller than 
that in SI, we expect the dynamics of the fluid evolution 
to be very similar to SO initially. As in the case of SI and 
S2, the outflow is expected to coUimate the magnetic field 
lines and generate magnetic shocks which may create tur- 
bulence in the torus. Turbulence can also be generated by 
the MRI, which operates on the orbital timescale of the 
torus independent of the field strength. We should then 
expect to see the stochastic accretion similar to the case 
in SI. Both a coUimated magnetic field and a massive, 
accretion torus surrounding a central black hole are es- 
sential ingredients for launching ultrarelativistic jets 56|. 
The black hole-torus system observed in our simulations 
provides a viable central engine for long-soft GRBs. 

The radial oscillation observed in model S2 gives rise 
to gravitational radiation. The oscillation period of 
^ 500Af corresponds to the gravitational wave frequency 
/ - l/[500M(l + z)] - 0.04(10''Mq/M)/(1 + z) Hz at 
redshift z. For a SMS with M > IO^'Mq, the signal 
is in the LISA frequency band. To estimate its ampli- 
tude, we apply the quadrupole formula h « 21- /Dl, 
where is the source's luminosity distance, f is the 
tracefree quadrupole moment and f ~ w^M^iskAi?^ 
2a;^Mdisk^-Rc- Here Rc ^ 30M is the characteristic ra- 
dius of the torus and A ^ 5M is the amplitude of the 
osciUation. Setting Mdisk 0.04M and uj = 2nf, we 
obtain 



4 X 10 



-23 



/ M 



V lO^Mp 



/48Gpc 



(35) 



where Dl = 48Gpc corresponds to redshift 2: = 5 in 
the concordance ACDM cosmology model with Hq = 
71km s-i Mpc-\ flM = 0.27 and = 0.73 We 
note that if the signal can be tracked for n cycles, where n 
is expected to be a few, the effective wave strength will be 
increased by a factor of ^/n. Such a gravitational wave 
signal may be detectable by LISA [sec [58] for LISA's 
sensitivity curve]. 

Our simulations are adiabatic and do not take into ac- 
count the heat loss due to neutrino cooling. To determine 
if this effect can be neglected during the phase in which 
the torus forms and evolves around the black hole, we es- 
timate the neutrino cooling timescale. We first compute 

the temperature in the disk from the specific thermal en- 

1/3 

ergy density eth = ^ - EcoM, where ecoid = ipo for our 
adopted F = 4/3 EOS. We find from our data that the 
typical values of po and eth in the disk at late times are 



Po 
eth/c^ 



6000 



/ M 



V 104 Afp 



g cm 



0.005 



(36) 
(37) 
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where we have restored the speed of hght c in the above 
equation. To calculate the temperature T, we adopt the 
expression of eth{po,T) in [5§|: 



eth 



2mpC^ 



3kT fl + 3Xn 



+ / 



(38) 



where k is the Boltzmann constant, a is radiation 
constant, nip is proton mass, and Xnuc is the mass 
fraction of free nucleons approximately given by [60j 
Xnuc « min[34.8p7(f/'*r{'/^exp(-0.61/Tn),l]. Here 
pio = /9o/10^"g cm-3 and Tn = T/10"if. The first 
term in Eq. psp is the specific thermal energy density 
of an ideal gas, and the second term is the contribution 
from thermal radiation. The quantity / is a temperature- 
dependent numerical factor depending on the species of 
ultrarelativistic particles that contribute to thermal ra- 
diation. When T > 2meC^/k 10^°/^, photons, ultra- 
relativistic electrons and positrons are present (assuming 
thermal neutrinos are negligible) and / — 11/4. On the 
other hand, when T <C lO^'^K, only photons contribute 
to thermal radiation and / = 1. Combining Eqs. (j36p - 



, we obtain 



0.0345(1 + 3Xnuc)Tg + 1.40/ 



M 



5 , (39) 



where Tg = T/\(fK. For M = IQ-^Mq, we find T w 
1.4 X IQ^K and, not surprisingly, eth is dominated by 
thermal photon radiation. At this density and tempera- 
ture, the torus is optically thin to neutrinos. The cooling 
rate Qi, is dominated by the pair neutrino process and the 
value is Qi, ~ lO^^erg cm~^ s~^ [6ll|. The neutrino cool- 
ing timescale is ^ Pa^th/Qv ^ 3 x lO^s ~ 5 x 10 ''A/, 
which is much longer than the timescale in our simu- 
lations (- 2000M). Even for M = lOOM©, we find 
Ti. - 90s - 2 X lO^Af > 2000M. The same conclu- 
sion (i.e. > 2000M) holds for all M > lOOMg. Hence 
neutrino cooling can be neglected in the torus evolution. 



V. SUMMARY AND CONCLUSION 

In this paper, we study the magnetorotational collapse 
of very massive stars by performing full GRMHD simu- 
lations in axisymmetry. We model the pre-coUapse star 
by an n = 3 polytrope uniformly rotating near the mass- 
shedding limit at the onset of radial collapse. We adopt 
an adiabatic F = 4/3 EOS for the fluid. We study three 
models, which we label SO, SI and S2. The three models 
differ by the strength of the initial magnetic field (see Ta- 
ble |T|. Model SO is unmagnctizcd [Ai — 0), whereas the 
ratio of the initial magnetic to kinetic energies [M /T) 
are 1% and 10% for models SI and S2, respectively. 

We find that these magnetic fields do not affect the ini- 
tial collapse significantly. An apparent horizon forms at 
time t K, 29000A'/. The black hole grows as the collapse 
proceeds, and settles down at a time ~ 150Af after the 



formation of the apparent horizon. For all three models 
we study, we find that the mass and spin parame- 
ter Jh/M'^ of the black hole are approximately 0.95Af 
and 0.7 respectively, where M is the initial mass of the 
star. These values roug hly agree with the semi-analytic 
estimates in [2^ [2^ l30l |. The remaining material forms a 
torus around the central black hole. Although the central 
black hole has settled down to quasi-stationary equilib- 
rium, the ambient torus continues to evolve as fluid from 
the outer layers of the star gradually reaches the central 
region. During this epoch, magnetic fields have substan- 
tial influence on the evolution of the torus. The infalling 
fluid particles have large angular momenta. They pile up 
near the black hole horizon, are heated by shocks and 
then get ejected along the surface of the torus, form- 
ing an unbound outflow. In the presence of magnetic 
flelds, the outflow bends the magnetic fleld lines near the 
boundary of the outflow, which ampliflcs the fleld and 
induces magnetic shocks. The extra magnetic pressure 
makes the outflow stronger than in the unmagnetized 
case. The outflow also causes the magnetic fields to col- 
limate along the black hole's rotation axis. For model 
SO, when the outflow leaves the central region, the torus 
settles down to equilibrium. For model SI, MHD turbu- 
lence generated by magnetic shocks and MRI in the disk 
causes stochastic accretion of material into the black hole. 
For model S2, when the outflow leaves, strong magnetic 
flelds in the torus create a magnetic wind, driving more 
material and magnetic fleld out of the torus. During this 
time, the torus acquires a quasiperiodic radial oscilla- 
tion. The wind subsides as the magnetic fleld inside the 
torus decreases. The radial oscillations of the torus in- 
duce episodic accretion of material into the central black 
hole. The oscillations also generate gravitational radia- 
tion, which might be detectable by LISA at redshift z ~ 5 
if the mass of the star satisfles M > IO^A/q. 

If the initial magnetic field strength is smaller than 
that in model SI, we expect the evolution to be similar 
to SI. In particular, the evolution in the collapse phase 
should remain unchanged. We also expect the outflow to 
collimate the magnetic fleld lines and generate magnetic 
shocks, which then leads to turbulence in the disk. Tur- 
bulence will be maintained as a result of the MRI. We 
thus expect stochastic accretion of the torus as in the 
case of SI. 

In typical cases, the final stage of the magnetorota- 
tional collapse consists of a central black hole surrounded 
by a coUimated magnetic field and a massive torus. These 
are the main ingredients for generating ultrarelativistic 
jets at large distance from the central source. The final 
system obtained in our simulations is thus capable of gen- 
erating a long-soft GRB. In principle, the collapse of a 
very massive star could result in the simultaneous detec- 
tion of gravitational waves and a GRB. The gravitational 
wave signal consists of an initial burst signal due to col- 
lapse, a black-hole ring-down signal, and a quasi-periodic 
signal due to the torus's oscillation if the magnetic field 
is strong. 
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A few issues warrant further study. The first is the 
EOS. Our r = 4/3 adiabatic EOS is a good approxi- 
mation only for very massive stars. But most of the ob- 
served long-soft GRBs are believed to be triggered by the 
magnetorotational core collapse of smaller-mass Pop I/II 
stars [sj]- The core mass of a Pop I/II star is less than 
2Mq. a r = 4/3 EOS describes the early phase of core 
collapse in such a star, when the pressure is dominated 
by relativistic degenerate electrons. But the EOS stiffens 
when the core density exceeds nuclear density and this 
happens before an apparent horizon forms. Also, a realis- 
tic EOS for this scenario must incorporate more detailed 
microphysics and neutrino transport. 

A second issue concerns a search for a more robust 
singularity- avoiding algorithm once a black hole forms. 
As mentioned in Sec. IIVI we are only able to evolve the 
system for ~ 200Af after the black hole formation with 
our current excision technique. However, the evolution 
timescale of the torus is > 2000M. While this evolution 
could be reliably tracked in the Cowling approximation, 
we are interested in more general scenarios. We plan to 
explore this issue in two ways. The first will be to search 
for better lapse and shift conditions that can suppress 
troublesome superluminal gauge modes. The other will 
be to identify a gauge that can drive the metric inside the 
horizon to a puncture-like solution, a technique which 
has been used with great success in binary black hole 
simulations (6^ . Simple experimentation with vacuum 
black holes and black holes immersed in hydrodynamic 
fluid suggest that there exist such gauge choices [63] . 

The third issue concerns our assumption of axisym- 
metry. Nonaxisymmetric instabilities such as bar and/or 



one-armed spiral instabilities may develop during the col- 
lapse, which could affect the subsequent dynamics (but 
see [6J] for a treatment of unmagnetized collapse in full 
3-1-1 post-Newtonian gravitation). Additionally, the 
MHD turbulence developed as a result of magnetic shocks 
and the MRI will be different. In particular, turbulence 
arises and persists more readily in 3 -I- 1 due to the lack 
of symmetry. More specifically, according to the axisym- 
metric anti-dynamo theorem [65| , sustained growth of the 
magnetic field energy is not possible through axisymmet- 
ric turbulence. However, a full 3-1-1 GRMHD simulation 
covering the required dynamic range for massive stellar 
collapse is computationally challenging and possibly be- 
yond the resources currently available. This is because 
the torus extends to a large distance away from the cen- 
tral black hole, requiring vast dynamic range, and the 
dynamical timescale of the torus is very long. Though 
simulations in full 3-1-1 dimensions will eventually be 
necessary to capture the full behavior of the collapse, the 
2-1-1 results presented here likely provide a reasonable 
first approximation. 
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